# ------------------------------------------------------------------------------
# Builds cost effectiveness variables
# From: stress_test_medicare repo <https://gitlab.com/labsysmed/zolab-projects/stress_test_medicare/-/blob/cabg_outcomes/code/03_analysis/02_prep-holdout/01_prep-holdout.R>
# Updates author: Cassidy Shubatt <cshubatt@gmail.com>
# To run: bash 11_cost_effectiveness.sh
# ------------------------------------------------------------------------------

# Seed ----------------------------------------------------------
set.seed(4815)

# Libraries ----------------------------------------------------------
library(here)
library(yaml)
library(data.table)
library(dplyr)
library(readr) # read_csv()
library(tidyr) # replace_na()
library(purrr) # walk()
library(glue) # glue()
library(feather) # read_feather()

u <- modules::use(here::here("lib", "util.R"))
temp <- here("code", "02_prep_and_summarize_cohort", "temp")

# Helpers ----------------------------------------------------------------------
get_daly <- function(d_int, stress, cath, int, int_days_to, const, ly_remaining) {
  early_int <- int & int_days_to <= 3

  # exp_lifeyears <- const$e_x

  p_benefit <- d_int *
    (stress | cath) *
    (0.82 * int + 0.18 * early_int)
  y_benefit <- ly_remaining * (const$p_fatal + (const$p_nonfatal * const$w_nonfatal))

  p_cost <- cath * const$r_cath_stroke
  y_cost <- const$w_stroke * ly_remaining

  return(p_benefit * y_benefit - p_cost * y_cost)
}

# Config -----------------------------------------------------------------------
message("Configuring...")
overnight_lab <- ""
paths <- read_yaml(here::here("lib", "filepaths.yml"))
const <- read_yaml(
  here("lib", "costeff_constants.yml")
)

# Load Data --------------------------------------------------------------------
message("Loading data...")
exp_lifeyears <- readRDS(paths$analysis$life_expectancy)
cohort <- readRDS(glue(paths$analysis$full_cohort)) %>%
  select(
    ptid, ed_enc_id, age_at_admit, stent_010_day, cabg_010_day, test_010_day,
    stress_010_day, cath_010_day
  ) %>%
  u$safe_left_join(exp_lifeyears)
days_to_int <- read_csv(paths$cohort$outcomes_p0) %>%
  mutate(
    days_to_stent = `days_to_stent-start_date_p0-start_date_p180`,
    days_to_cabg = `days_to_cabg-start_date_p0-start_date_p180`
  ) %>%
  select(ed_enc_id, days_to_stent, days_to_cabg)

# Cost Effectiveness -----------------------------------------------------------
message("Calculating cost effectiveness...")
visits <- cohort %>%
  u$safe_left_join(days_to_int) %>%
  setDT()

visits[, cost := (const$cost_effectiveness$c_stress * stress_010_day) +
  (const$cost_effectiveness$c_cath * cath_010_day) +
  (const$cost_effectiveness$c_stent * stent_010_day) +
  (const$cost_effectiveness$c_cabg * cabg_010_day)]

walk(
  const$cost_effectiveness$d_int,
  ~ visits[, paste0("daly_", 1000 * .x) := get_daly(
    .x,
    stress_010_day,
    cath_010_day,
    stent_010_day | cabg_010_day,
    pmin(days_to_stent, days_to_cabg, Inf, na.rm = TRUE),
    const$cost_effectiveness, ly_remaining
  )]
)

# Save -------------------------------------------------------------------------
save_fp <- paths$analysis$daly_cost
message("Saving daly df to ", save_fp, "...")
write_rds(visits, save_fp)

message("Done.")
